//
// File: calc_TimeFreq_Lomb.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 21-Mar-2025 16:18:04
//

// Include Files
#include "calc_TimeFreq_Lomb.h"
#include "blockedSummation.h"
#include "calcAreas.h"
#include "calc_TimeFreq_Lomb_types.h"
#include "detrend.h"
#include "diff.h"
#include "linspace.h"
#include "lomb2.h"
#include "mean.h"
#include "nullAssignment.h"
#include "qrsolve.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include "rt_defines.h"
#include <cmath>

// Function Declarations
static double rt_atan2d_snf(double u0, double u1);

// Function Definitions
//
// Arguments    : double u0
//                double u1
// Return Type  : double
//
static double rt_atan2d_snf(double u0, double u1)
{
  double y;
  if (std::isnan(u0) || std::isnan(u1)) {
    y = rtNaN;
  } else if (std::isinf(u0) && std::isinf(u1)) {
    int b_u0;
    int b_u1;
    if (u0 > 0.0) {
      b_u0 = 1;
    } else {
      b_u0 = -1;
    }
    if (u1 > 0.0) {
      b_u1 = 1;
    } else {
      b_u1 = -1;
    }
    y = std::atan2(static_cast<double>(b_u0), static_cast<double>(b_u1));
  } else if (u1 == 0.0) {
    if (u0 > 0.0) {
      y = RT_PI / 2.0;
    } else if (u0 < 0.0) {
      y = -(RT_PI / 2.0);
    } else {
      y = 0.0;
    }
  } else {
    y = std::atan2(u0, u1);
  }
  return y;
}

//
// Arguments    : const coder::array<double, 1U> &t
//                const coder::array<double, 1U> &y
//                double fs
//                double nfft
//                double winSize
//                double overlap
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct0_T *lomb_out
// Return Type  : void
//
void calc_TimeFreq_Lomb(const coder::array<double, 1U> &t,
                        const coder::array<double, 1U> &y, double fs,
                        double nfft, double winSize, double overlap,
                        const double VLF[2], const double LF[2],
                        const double HF[2], struct0_T *lomb_out)
{
  coder::array<double, 2U> W;
  coder::array<double, 2U> idx;
  coder::array<double, 2U> r;
  coder::array<double, 1U> c_y;
  coder::array<double, 1U> iStart;
  coder::array<double, 1U> wtmt;
  coder::array<double, 1U> y2;
  coder::array<double, 1U> z;
  coder::array<int, 2U> b_y;
  coder::array<int, 1U> c_i;
  coder::array<boolean_T, 1U> x;
  struct2_T expl_temp;
  double b_t;
  double w;
  int b_i;
  int b_iStart;
  int cnt;
  int i;
  int i1;
  int iNext;
  int loop_ub;
  int nSeg;
  w = fs / 2.0;
  // calLomb - Calculates PSD using windowed Lomb-Scargle method.
  //
  // Inputs:
  // Outputs:
  // get limits of windows
  if (t[t.size(0) - 1] >= winSize) {
    boolean_T flag;
    iStart.set_size(t.size(0));
    loop_ub = t.size(0);
    for (i = 0; i < loop_ub; i++) {
      iStart[i] = t[i];
    }
    // windowArray - Seperates the input array s into smaller segments/windows
    // with lengths of winSize. winSize can be either in units of time or
    // samples. If winSize is given in units of time the first column of the
    // input array must containg time and method should be set to a value of 1.
    //
    // INPUTS:
    //    s: input array OR length of array. If outputing segs or using the time
    //        methods, the 1st dim must contain time values.
    //    winSize: size of segments to return
    //    method: method of segmenting array (0=by samples, 1=by time)
    //    overlap: amount that segments/windows overlap.
    //    flagCenters: decides whether to output center of windows instead of
    //    beginning
    //
    // OUTPUTS:
    //    idx: 2dim array that contains the start and end of each segment/window
    //    centers: array containg locations of center of windows (useful when
    //    using the time method)
    //
    // EXAMPLES:
    //    1. Seperate input using window size of 100 seconds and overlap
    //    of 50 seconds.
    //    idx=slidingWindow(input,100,50,1);
    //
    //    2. Seperate input using window size of 128 samples and overlap
    //    of 64 samples.
    //    idx=slidingWindow(input,128,64,0);
    //
    //    3. Seperate input using window size of 128 samples and overlap
    //    of 64 samples. Also outputs a cell array containing the actual
    //    segments [idx, segs]=slidingWindow(input,128,64,0,1);
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
    //
    //  This file is part of HRVAS
    //
    //  HRVAS is free software: you can redistribute it and/or modify
    //  it under the terms of the GNU General Public License as published by
    //  the Free Software Foundation, either version 3 of the License, or
    //  (at your option) any later version.
    //
    //  HRVAS is distributed in the hope that it will be useful,
    //  but WITHOUT ANY WARRANTY; without even the implied warranty of
    //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    //  GNU General Public License for more details.
    //
    //  You should have received a copy of the GNU General Public License
    //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // Check inputs
    // includeSeg=0;
    // check inputs
    //  if (m==1) && (n==1) && (method==1)
    //      error('slidingWindow: Invalid input array for this windowing
    //      method.') return;
    //  end
    //
    //  if (m==1) && (n==1)
    //      N=s;
    //  %     if includeSeg
    //  %         warning('SlidingWindow: Cannot output segs if input is
    //  scalar'); %     end
    //      includeSeg=false;
    //  else
    //      N=length(s);
    //  end
    // segment/window data based on time
    nSeg = static_cast<int>(std::floor(static_cast<double>(t.size(0)) / 2.0));
    // don't know how many segments there will be...so assume no more than N/2
    cnt = -1;
    b_i = 0;
    b_iStart = 1;
    iNext = 1;
    flag = true;
    idx.set_size(nSeg, 2);
    loop_ub = nSeg << 1;
    for (i = 0; i < loop_ub; i++) {
      idx[i] = 0.0;
    }
    // preallocate idx array
    loop_ub = t.size(0);
    for (i = 0; i < loop_ub; i++) {
      iStart[i] = t[i] - t[0];
    }
    // offset all to make 1st time point = 0
    if (overlap == 0.0) {
      boolean_T exitg1;
      // if overlap is zero use this method to get win limits...faster
      // time
      b_iStart = iStart.size(0) - 1;
      for (i = 0; i <= b_iStart; i++) {
        iStart[i] = iStart[i] / 5.0;
      }
      iStart.set_size(b_iStart + 1);
      b_iStart++;
      for (int k{0}; k < b_iStart; k++) {
        iStart[k] = std::floor(iStart[k]);
      }
      //
      coder::diff(iStart, c_y);
      x.set_size(c_y.size(0));
      loop_ub = c_y.size(0);
      for (i = 0; i < loop_ub; i++) {
        x[i] = (c_y[i] == 1.0);
      }
      b_iStart = x.size(0);
      nSeg = 0;
      c_i.set_size(x.size(0));
      cnt = 0;
      exitg1 = false;
      while ((!exitg1) && (cnt <= b_iStart - 1)) {
        if (x[cnt]) {
          nSeg++;
          c_i[nSeg - 1] = cnt + 1;
          if (nSeg >= b_iStart) {
            exitg1 = true;
          } else {
            cnt++;
          }
        } else {
          cnt++;
        }
      }
      if (x.size(0) == 1) {
        if (nSeg == 0) {
          c_i.set_size(0);
        }
      } else {
        if (nSeg < 1) {
          nSeg = 0;
        }
        c_i.set_size(nSeg);
      }
      lomb_out->global_PSD.set_size(c_i.size(0));
      loop_ub = c_i.size(0);
      for (i = 0; i < loop_ub; i++) {
        lomb_out->global_PSD[i] = c_i[i];
      }
      // locations of the end of all windows
      iStart.set_size(lomb_out->global_PSD.size(0) + 1);
      iStart[0] = 1.0;
      loop_ub = lomb_out->global_PSD.size(0);
      for (i = 0; i < loop_ub; i++) {
        iStart[i + 1] = lomb_out->global_PSD[i] + 1.0;
      }
      if (iStart.size(0) - 1 < 1) {
        loop_ub = 0;
      } else {
        loop_ub = iStart.size(0) - 1;
      }
      idx.set_size(loop_ub, 2);
      for (i = 0; i < loop_ub; i++) {
        idx[i] = iStart[i];
      }
      loop_ub = lomb_out->global_PSD.size(0);
      for (i = 0; i < loop_ub; i++) {
        idx[i + idx.size(0)] = lomb_out->global_PSD[i];
      }
    } else {
      // overlap > 0 ... use loops
      while (b_i <= t.size(0) - 1) {
        b_i++;
        b_t = iStart[b_i - 1] - iStart[b_iStart - 1];
        // add time values
        if (flag && (winSize - b_t <= overlap)) {
          iNext = b_i;
          flag = false;
        } else if (b_t >= winSize) {
          cnt++;
          idx[cnt] = b_iStart;
          idx[cnt + idx.size(0)] = b_i - 1;
          b_i = iNext - 1;
          b_iStart = iNext;
          flag = true;
        }
      }
      b_y.set_size(1, (idx.size(0) - cnt) - 1);
      loop_ub = (idx.size(0) - cnt) - 2;
      for (i = 0; i <= loop_ub; i++) {
        b_y[i] = (cnt + i) + 2;
      }
      coder::internal::nullAssignment(idx, b_y);
      // remove extra zeros created during preallocation
    }
    //  if includeSeg
    //      segs=cell(size(idx,1));
    //      for i=1:size(idx,1)
    //          segs{i}=s(idx(i,1):idx(i,2));%code here for outputing
    //      end
    //  end
  } else {
    idx.set_size(1, 2);
    idx[0] = 1.0;
    idx[1] = t.size(0);
  }
  b_t = std::round(winSize / 2.0);
  loop_ub = idx.size(0);
  lomb_out->T.set_size(idx.size(0));
  for (i = 0; i < loop_ub; i++) {
    lomb_out->T[i] = t[static_cast<int>(idx[i]) - 1] + b_t;
  }
  // estimate the center of the windows for plotting
  // preallocate memory
  // number of PSD/windows
  lomb_out->psd.set_size(static_cast<int>(nfft), idx.size(0));
  coder::linspace(w - w / nfft, nfft, r);
  lomb_out->f.set_size(r.size(1));
  loop_ub = r.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->f[i] = r[i];
  }
  i = idx.size(0);
  for (b_i = 0; b_i < i; b_i++) {
    int i2;
    int i3;
    int i4;
    // Prepare y2 and t2
    b_t = idx[b_i];
    w = idx[b_i + idx.size(0)];
    if (b_t > w) {
      i1 = 0;
      i2 = -1;
      i3 = 0;
      i4 = 0;
    } else {
      i1 = static_cast<int>(b_t) - 1;
      i2 = static_cast<int>(w) - 1;
      i3 = static_cast<int>(b_t) - 1;
      i4 = static_cast<int>(w);
    }
    // remove linear trend
    loop_ub = i2 - i1;
    y2.set_size(loop_ub + 1);
    for (b_iStart = 0; b_iStart <= loop_ub; b_iStart++) {
      y2[b_iStart] = y[i1 + b_iStart];
    }
    if (loop_ub < 0) {
      nSeg = 0;
    } else {
      nSeg = loop_ub + 1;
    }
    b_y.set_size(1, nSeg);
    if (nSeg > 0) {
      b_y[0] = 0;
      cnt = 0;
      for (int k{2}; k <= nSeg; k++) {
        cnt++;
        b_y[k - 1] = cnt;
      }
    }
    iStart.set_size(b_y.size(1));
    nSeg = b_y.size(1);
    for (b_iStart = 0; b_iStart < nSeg; b_iStart++) {
      iStart[b_iStart] = b_y[b_iStart];
    }
    if (loop_ub + 1 != 0) {
      if (loop_ub + 1 == 1) {
        y2.set_size(loop_ub + 1);
        for (i2 = 0; i2 <= loop_ub; i2++) {
          y2[i2] = y[i1 + i2] * 0.0;
        }
      } else {
        double p[2];
        cnt = iStart.size(0) - 1;
        lomb_out->global_PSD.set_size(iStart.size(0));
        for (iNext = 0; iNext <= cnt; iNext++) {
          lomb_out->global_PSD[iNext] = iStart[iNext] - iStart[0];
        }
        nSeg = lomb_out->global_PSD.size(0);
        for (b_iStart = 0; b_iStart < nSeg; b_iStart++) {
          lomb_out->global_PSD[b_iStart] =
              lomb_out->global_PSD[b_iStart] / iStart[iStart.size(0) - 1];
        }
        nSeg = lomb_out->global_PSD.size(0);
        for (b_iStart = 0; b_iStart < nSeg; b_iStart++) {
          b_t = lomb_out->global_PSD[b_iStart];
          lomb_out->global_PSD[b_iStart] = std::fmax(b_t, 0.0);
        }
        W.set_size(iStart.size(0), 2);
        for (iNext = 0; iNext <= cnt; iNext++) {
          W[iNext] = lomb_out->global_PSD[iNext];
        }
        for (iNext = 0; iNext <= cnt; iNext++) {
          W[iNext + W.size(0)] = 1.0;
        }
        c_y.set_size(loop_ub + 1);
        for (b_iStart = 0; b_iStart <= loop_ub; b_iStart++) {
          c_y[b_iStart] = y[i1 + b_iStart];
        }
        coder::internal::qrsolve(W, c_y, p, &nSeg);
        nSeg = W.size(0);
        y2.set_size(W.size(0));
        for (iNext = 0; iNext < nSeg; iNext++) {
          y2[iNext] = W[iNext] * p[0] + W[W.size(0) + iNext] * p[1];
        }
        if (loop_ub + 1 == y2.size(0)) {
          y2.set_size(loop_ub + 1);
          for (i2 = 0; i2 <= loop_ub; i2++) {
            y2[i2] = y[i1 + i2] - y2[i2];
          }
        } else {
          binary_expand_op(y2, y, i1, i2);
        }
      }
    }
    b_t = coder::blockedSummation(y2, y2.size(0)) /
          static_cast<double>(y2.size(0));
    loop_ub = y2.size(0);
    for (i1 = 0; i1 < loop_ub; i1++) {
      y2[i1] = y2[i1] - b_t;
    }
    // remove mean
    // Calculate un-normalized lomb PSD
    //
    //   [Pn] = lomb(t, y, f)
    //
    //   Uses Lomb's method to compute normalized
    //   periodogram values "Pn" as a function of
    //   supplied vector of frequencies "f" for
    //   input vectors "t" (time) and "y" (observations).
    //   Also returned is probability "Prob" of same
    //   length as Pn (and f) that the null hypothesis
    //   is valid. If f is not supplied it assumes
    //   f =  [1/1024 : 1/1024 : 0.5/min(diff(t))];
    //   x and y must be the same length.
    //  See also:
    //
    //  [1] N.R. Lomb, ``Least-squares frequency analysis of
    //  unequally spaced data,'' Astrophysical and Space Science,
    //  (39) pp. 447--462, 1976.   ... and
    //
    //  [2] J.~D. Scargle, ``Studies in astronomical time series analysis.
    //  II. Statistical aspects of spectral analysis of unevenly spaced data,''
    //  Astrophysical Journal, vol. 263, pp. 835--853, 1982.
    //
    //  [3] T. Thong, "Lomb-Welch Periodogram for Non-uniform Sampling",
    //  Proceedings for the 26th anual international conference of the IEEE
    //  EMBS, Sept 1-5, 2004. check inputs subtract mean, compute variance,
    //  initialize Pn
    b_t = coder::blockedSummation(y2, y2.size(0)) /
          static_cast<double>(y2.size(0));
    z.set_size(y2.size(0));
    loop_ub = y2.size(0);
    for (i1 = 0; i1 < loop_ub; i1++) {
      z[i1] = y2[i1] - b_t;
    }
    c_y.set_size(lomb_out->f.size(0));
    loop_ub = lomb_out->f.size(0);
    for (i1 = 0; i1 < loop_ub; i1++) {
      c_y[i1] = 0.0;
    }
    // 	now do main loop for all frequencies
    i1 = lomb_out->f.size(0);
    for (iNext = 0; iNext < i1; iNext++) {
      w = 6.2831853071795862 * lomb_out->f[iNext];
      if (w > 0.0) {
        double a;
        a = 2.0 * w;
        loop_ub = i4 - i3;
        lomb_out->global_PSD.set_size(loop_ub);
        for (i2 = 0; i2 < loop_ub; i2++) {
          lomb_out->global_PSD[i2] = a * t[i3 + i2];
        }
        iStart.set_size(lomb_out->global_PSD.size(0));
        nSeg = lomb_out->global_PSD.size(0);
        for (i2 = 0; i2 < nSeg; i2++) {
          iStart[i2] = lomb_out->global_PSD[i2];
        }
        b_iStart = lomb_out->global_PSD.size(0);
        for (int k{0}; k < b_iStart; k++) {
          iStart[k] = std::sin(iStart[k]);
        }
        b_iStart = lomb_out->global_PSD.size(0);
        for (int k{0}; k < b_iStart; k++) {
          lomb_out->global_PSD[k] = std::cos(lomb_out->global_PSD[k]);
        }
        b_t = rt_atan2d_snf(
                  coder::blockedSummation(iStart, iStart.size(0)),
                  coder::blockedSummation(lomb_out->global_PSD,
                                          lomb_out->global_PSD.size(0))) /
              2.0 / w;
        wtmt.set_size(loop_ub);
        for (i2 = 0; i2 < loop_ub; i2++) {
          wtmt[i2] = w * (t[i3 + i2] - b_t);
        }
        iStart.set_size(wtmt.size(0));
        loop_ub = wtmt.size(0);
        for (i2 = 0; i2 < loop_ub; i2++) {
          iStart[i2] = wtmt[i2];
        }
        b_iStart = wtmt.size(0);
        for (int k{0}; k < b_iStart; k++) {
          iStart[k] = std::cos(iStart[k]);
        }
        lomb_out->global_PSD.set_size(wtmt.size(0));
        loop_ub = wtmt.size(0);
        for (i2 = 0; i2 < loop_ub; i2++) {
          lomb_out->global_PSD[i2] = wtmt[i2];
        }
        b_iStart = wtmt.size(0);
        for (int k{0}; k < b_iStart; k++) {
          lomb_out->global_PSD[k] = std::sin(lomb_out->global_PSD[k]);
        }
        loop_ub = z.size(0);
        if (z.size(0) == iStart.size(0)) {
          iStart.set_size(z.size(0));
          for (i2 = 0; i2 < loop_ub; i2++) {
            iStart[i2] = z[i2] * iStart[i2];
          }
        } else {
          times(iStart, z);
        }
        a = coder::blockedSummation(iStart, iStart.size(0));
        loop_ub = z.size(0);
        if (z.size(0) == lomb_out->global_PSD.size(0)) {
          lomb_out->global_PSD.set_size(z.size(0));
          for (i2 = 0; i2 < loop_ub; i2++) {
            lomb_out->global_PSD[i2] = z[i2] * lomb_out->global_PSD[i2];
          }
        } else {
          binary_expand_op(lomb_out, z);
        }
        w = coder::blockedSummation(lomb_out->global_PSD,
                                    lomb_out->global_PSD.size(0));
        iStart.set_size(wtmt.size(0));
        loop_ub = wtmt.size(0);
        for (i2 = 0; i2 < loop_ub; i2++) {
          iStart[i2] = wtmt[i2];
        }
        b_iStart = wtmt.size(0);
        for (int k{0}; k < b_iStart; k++) {
          iStart[k] = std::cos(iStart[k]);
        }
        b_iStart = wtmt.size(0);
        for (int k{0}; k < b_iStart; k++) {
          wtmt[k] = std::sin(wtmt[k]);
        }
        loop_ub = iStart.size(0);
        for (i2 = 0; i2 < loop_ub; i2++) {
          b_t = iStart[i2];
          iStart[i2] = b_t * b_t;
        }
        lomb_out->global_PSD.set_size(wtmt.size(0));
        loop_ub = wtmt.size(0);
        for (i2 = 0; i2 < loop_ub; i2++) {
          b_t = wtmt[i2];
          lomb_out->global_PSD[i2] = b_t * b_t;
        }
        c_y[iNext] = a * a / coder::blockedSummation(iStart, iStart.size(0)) +
                     w * w /
                         coder::blockedSummation(lomb_out->global_PSD,
                                                 lomb_out->global_PSD.size(0));
      } else {
        double a;
        loop_ub = z.size(0);
        nSeg = i4 - i3;
        if (z.size(0) == nSeg) {
          iStart.set_size(z.size(0));
          for (i2 = 0; i2 < loop_ub; i2++) {
            iStart[i2] = z[i2] * t[i3 + i2];
          }
        } else {
          binary_expand_op(iStart, z, t, i3, i4 - 1);
        }
        a = coder::blockedSummation(iStart, iStart.size(0));
        iStart.set_size(nSeg);
        for (i2 = 0; i2 < nSeg; i2++) {
          b_t = t[i3 + i2];
          iStart[i2] = b_t * b_t;
        }
        c_y[iNext] = a * a / coder::blockedSummation(iStart, iStart.size(0));
      }
    }
    //  return denormalized spectrum (see T. Thong)
    loop_ub = c_y.size(0);
    for (i1 = 0; i1 < loop_ub; i1++) {
      lomb_out->psd[i1 + lomb_out->psd.size(0) * b_i] =
          c_y[i1] / static_cast<double>(y2.size(0));
    }
  }
  //  calcAreas - Calulates areas/energy under the PSD curve within the freq
  //  bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
  //  percentage, and normalized units. Also returns LF/HF ratio.
  //
  //  Inputs:
  //    PSD: PSD vector
  //    F: Freq vector
  //    VLF, LF, HF: array containing VLF, LF, and HF freq limits
  //    flagVLF: flag to decide whether to calculate VLF hrv
  //  Output:
  //
  //  Usage:
  //
  //
  //  Ref: Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
  lomb_out->hrv.aVLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.aVLF[i] = 0.0;
  }
  lomb_out->hrv.aLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.aLF[i] = 0.0;
  }
  lomb_out->hrv.aHF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.aHF[i] = 0.0;
  }
  lomb_out->hrv.aTotal.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.aTotal[i] = 0.0;
  }
  lomb_out->hrv.pVLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.pVLF[i] = 0.0;
  }
  lomb_out->hrv.pLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.pLF[i] = 0.0;
  }
  lomb_out->hrv.pHF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.pHF[i] = 0.0;
  }
  lomb_out->hrv.nLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.nLF[i] = 0.0;
  }
  lomb_out->hrv.nHF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.nHF[i] = 0.0;
  }
  lomb_out->hrv.LFHF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.LFHF[i] = 0.0;
  }
  lomb_out->hrv.peakVLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.peakVLF[i] = 0.0;
  }
  lomb_out->hrv.peakLF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.peakLF[i] = 0.0;
  }
  lomb_out->hrv.peakHF.set_size(lomb_out->psd.size(1));
  loop_ub = lomb_out->psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->hrv.peakHF[i] = 0.0;
  }
  i = lomb_out->psd.size(1);
  for (nSeg = 0; nSeg < i; nSeg++) {
    loop_ub = lomb_out->psd.size(0);
    c_y.set_size(lomb_out->psd.size(0));
    for (i1 = 0; i1 < loop_ub; i1++) {
      c_y[i1] = lomb_out->psd[i1 + lomb_out->psd.size(0) * nSeg];
    }
    calcAreas(lomb_out->f, c_y, VLF, LF, HF, &expl_temp);
    // create output structure
    lomb_out->hrv.aVLF[nSeg] = expl_temp.aVLF;
    lomb_out->hrv.aLF[nSeg] = expl_temp.aLF;
    lomb_out->hrv.aHF[nSeg] = expl_temp.aHF;
    lomb_out->hrv.aTotal[nSeg] = expl_temp.aTotal;
    lomb_out->hrv.pVLF[nSeg] = expl_temp.pVLF;
    lomb_out->hrv.pLF[nSeg] = expl_temp.pLF;
    lomb_out->hrv.pHF[nSeg] = expl_temp.pHF;
    lomb_out->hrv.nLF[nSeg] = expl_temp.nLF;
    lomb_out->hrv.nHF[nSeg] = expl_temp.nHF;
    lomb_out->hrv.LFHF[nSeg] = expl_temp.LFHF;
    lomb_out->hrv.peakVLF[nSeg] = expl_temp.peakVLF;
    lomb_out->hrv.peakLF[nSeg] = expl_temp.peakLF;
    lomb_out->hrv.peakHF[nSeg] = expl_temp.peakHF;
  }
  x.set_size(lomb_out->hrv.LFHF.size(0));
  loop_ub = lomb_out->hrv.LFHF.size(0);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (lomb_out->hrv.LFHF[i] > 1.0);
  }
  nSeg = x.size(0);
  if (x.size(0) == 0) {
    b_iStart = 0;
  } else {
    b_iStart = x[0];
    for (int k{2}; k <= nSeg; k++) {
      b_iStart += x[k - 1];
    }
  }
  x.set_size(lomb_out->hrv.LFHF.size(0));
  loop_ub = lomb_out->hrv.LFHF.size(0);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (lomb_out->hrv.LFHF[i] <= 1.0);
  }
  nSeg = x.size(0);
  if (x.size(0) == 0) {
    cnt = 0;
  } else {
    cnt = x[0];
    for (int k{2}; k <= nSeg; k++) {
      cnt += x[k - 1];
    }
  }
  // global psd
  coder::mean(lomb_out->psd, lomb_out->global_PSD);
  calcAreas(lomb_out->f, lomb_out->global_PSD, VLF, LF, HF,
            &lomb_out->global_hrv);
  lomb_out->hrv.rLFHF =
      static_cast<double>(b_iStart) / static_cast<double>(cnt);
  lomb_out->global_f.set_size(lomb_out->f.size(0));
  loop_ub = lomb_out->f.size(0);
  for (i = 0; i < loop_ub; i++) {
    lomb_out->global_f[i] = lomb_out->f[i];
  }
}

//
// File trailer for calc_TimeFreq_Lomb.cpp
//
// [EOF]
//
